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We show that generic, slow dynamics can occur in the contact process on complex networks 
with a tree-like structure and a superimposed weight pattern, in the absence of additional (non- 
topological) sources of quenched disorder. The slow dynamics is induced by rare-region effects 
occurring on correlated subspaces of vertices connected by large weight edges, and manifests in the 
form of a smeared phase transition. We conjecture that more sophisticated network motifs could be 
able to induce Griffiths phases, as a consequence of purely topological disorder. 



PACS numbers: 89. 75. He. 05.70.Ln, 89.75.Fb 



I. INTRODUCTION 



The science of complex networks has witnessed a veri- 
table explosion of activity in the last decade, having be- 
come a powerful tool to represent, analyze and under- 
stand a myriad of natural and man-made systems, char- 
acterized by heterogeneous topological structures [l|, Q ■ 
In parallel with studies focused on the topology and func- 
tion of these systems, a large research effort has been de- 
voted to the analysis of non-equilibrium dynamical pro- 
cesses running on top of complex networks and on the 
effects that those complex substrates can have on their 
temporal behavior and also on possible phase transitions 
[H, Q . Such effects are particularly important in networks 
with a heterogeneous contact pattern, as in scalc-frcc 
(SF) networks @, in which the degree distribution, de- 
fined as the probability that an element (vertex or node) 
is connected to k others (has degree k), exhibits long tails 
as P{k) ~ fc~ 7 , with a certain degree exponent usually 
in the range 2 < 7 < 3 [10. 

One of the simplest non-equilibrium processes studied 
on networks is the contact process (CP) [1, 0- I n this 
model, vertices can be either occupied or empty. Empty 
vertices become occupied on contact with an occupied 
neighbor of degree hi with a rate A/fc^. On the other 
hand, occupied vertices become empty with a unitary 
rate. On a regular lattice the CP experiences a non- 
equilibrium phase transition at a critical point A^, sep- 
arating an absorbing phase from an active one [a-llCf. 
whose order parameter is the density of occupied sites 
p in the steady state. Thus, for A < A c , an absorb- 
ing phase with p = is observed, while for A > A c the 
system reaches an active phase, with p > in the ther- 
modynamic limit. Through a systematic analysis rely- 
ing on numerical simulations [Ill4l4 1 and theoretical ap- 
proaches based on the heterogeneous mean-field (HMF) 
theory [ll|, [T^, [l5[ a picture has emerged on the behav- 
ior of the CP on complex networks that emphasizes the 



strong effects of the network heterogeneity In particular, 
it has been shown that its absorbing pha se transition ex- 
hibits a nontrivial finite-size scaling jig , depending not 
only on the number of vertices N, but also on the degree 
fluctuations of the network, measured by the second mo- 
ment of the degree distribution (A: 2 ) = ^2 k k 2 P(k). This 
dependence induces very strong corrections to scaling in 
SF networks, which, if properly taken into account, can 
actually be observed in numerical simulations (l3l . ITij . 

The study of the CP, as well as other processes, both 
in and out of equilibrium 0, has shown the impor- 
tant effects that the disordered, heterogeneous topologi- 
cal structure in a network can impose on the dynamical 
systems defined on top of it. The effects of network dis- 
order have recently been extended one step further, in a 
series of papers [17Hl9{ dealing with the possibility of ob- 
serving Griffiths phases (GP) and rare-region (RR) phe- 
nomena 0, [2l[ at the phase transition of the CP on com- 
plex networks 1 . In case of regular lattices it is well known 
that quenched disorder 2 can strongly alter the behavior 
of a phase transition, imposing an anomalously slow re- 
laxation [23| ■ This is the consequence of inhomogeneities 
that can create, in the thermodynamic limit, RRs of char- 
acteristic size /, with probability Prr(1) oc exp(— cl), in 
which the system can stay active for exponentially long 
times t{1) cx cxp(W) below the critical point A c . In sys- 
tems with spatially uncorrelated disorder, a sharp phase 
transition is preserved, although it occurs at a critical 
point A c larger than the corresponding to the clean sys- 
tem \®. In the region A° < A < A c , a GP develops 
with a slow algebraic density decay p cx t~ c l h , c/b be- 



1 See also Ref. |22j for an investigation of GPs in the context of 
the superconductor-insulator transition on SF networks. 

2 When the disorder is annealed it acts as a noise and it is hence 
irrelevant for the universal scaling behavior of the phase transi- 
tion. 
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ing a non-universal exponent, which can be understood 
by non-perturbative methods j2^ - |28j . However, if the 
inhomogencitics arc correlated in a subspace with a di- 
verging diameter in dimensions d~ < cIrr < d, where d is 
the dimension of the system and d~ is the lower critical 
dimension, the transition becomes completely smeared. 
The smearing is due to fact that RRs can become infi- 
nite (in the thermodynamic limit) and undergo indepen- 
dent phase transitions, ordering at different values of A 
[2l| . This situation prevents the development of global 
order, in such a way that the clean sharp transition to 
the absorbing state in destroyed, and the associated sin- 
gularities become rounded. In this case, above the clean 
critical point, the density remains finite in the long time 
limit and approaches this value in a power-law form 



p(t) - p(oo) ~ t-< 



(1) 



where the exponent C is a ls° non-universal. Close to 
A°, however, the initial decay of the density follows an 
stretched exponential form 



hi[p(i)] cx & /{d ^ +z) , 



(2) 



where z is the dynamical exponent and d r is the dimen- 
sion of the uncorrclated subspace, as can be seen by ap- 
plying optimal fluctuation theory arguments plj . 

In the case of the CP on complex networks, it has re- 
cently been shown that an intrinsic quenched disorder 3 
defined by a varying control parameter Ai, can induce 
GPs and other RR effects on random Erdos-Renyi net- 
works [3(| below (and at) the percolation threshold (TtI — 
Il9j . The authors in [13, EH also considered the issue of 
the effects of pure topological disorder for a CP with con- 
stant A. By means of theoretical arguments and numeri- 
cal simulations, they conjectured that RR effects should 
be relevant only on networks with a finite topological di- 
mension D, defined in terms of the number of vertices at 
a topological distance £ from a given source, N(t) ~ £ D . 
Therefore, for random networks, in which the small- world 
property [3l[ holds, namely, where N(£) ~ exp(£) and 
D — > oo, RR effects and GPs should be absent. 

In this paper we reexamine the issue of rare-region ef- 
fects and the associated slow dynamics in the CP on SF 
small-world networks with additional topological disor- 
der beyond random connectivity. In particular, we con- 
sider the effect of two topological features which have 
been shown to induce a strong slowing down in dynami- 
cal processes, namely a tree-like structured network [Hj] 
and a weight pattern superimposed on it [33| . By means 
of extensive simulations, we provide numerical evidence 
that, although GPs cannot actually be observed, the 
loop-less tree-like structure together with an heteroge- 
neous weight pattern induces a smearing of the CP crit- 
ical phase transition and slow, power-law dynamics for 



finite times. This can be understood via RR effects trig- 
gered by topological correlations in the weights. Our 
results open the path to consider more complex topolog- 
ical motifs, such as a large clustering [3l| or a hierarchical 
community structure [341 ] . as generators of purely topo- 
logical Griffiths phases in non-equilibrium phase transi- 
tion on complex networks. 



II. THE CONTACT PROCESS ON WEIGHTED 
TREES 

Both a tree-like structure and the presence of a weight 
distribution are by themselves capable of slowing down 
the temporal evolution of dynamical processes on ran- 
dom networks [H, [H, HU . The induced slowing down is 
due to the generation of topological traps that capture 
dynamics and prevent the fast and wide range spread- 
ing naturally expected in small-world networks. In tree 
networks, traps are created by the lack of loops, which 
implies that there is a single path between any two ver- 
tices: Once activity is deep in the leaves of a subtree, it 
cannot reach other sections of the network until it first 
finds the exit from that subtree [32] • In weighted net- 
works, on the other hand, activity can become trapped 
in correlated sets of adjacent vertices, joined by edges 
carrying a large weight, which prevent the exploration of 
other regions of the network [33|. We expect that the 
combination of those slowing down ingredients can lead 
to the emergence of RR phenomena. 



A. Network models 

We consider the CP on weighted tree networks con- 
structed following the Barabasi- Albert (BA) model 
The choice of this model is motivated by the fact that it 
allows to construct tree structures in a very simple way, 
in contrast with other standard network generation mod- 
els, e.g Ref. [Hj]. The BA is a growing network model in 
which, at each time step s, a new vertex with m edges is 
added to the network and connected to an existing vertex 
s' of degree k s > with probability n s _> s / = k a ' /J2s"<s ^ s " • 
This process is iterated until reaching the desired network 
size N. The resulting network has a SF degree distribu- 
tion P(k) ~ 2m 2 fc~ 3 ; additionally, fixing m = 1 leads to 
a strictly tree (loop-less) topology. 

Binary (non-weighted) BA networks can be trans- 
formed into weighted ones by assigning to every edge 
connecting vertices i and j a symmetric weight cjy. We 
have considered two different schemes for weight assign- 
ment based on the topology of the network, i.e. weights 
are not random but depend on properties of the vertices 
they connect. The two strategies are defined as follows: 

(i) Weighted BA tree I (WBAT-I): Multiplicative 
weights depending on the degree of the adjacent vertices, 
namely 



Temporal GPs have been considered in Ref. [29 



Wij = LU {kikj)~ 



(3) 
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where uq is an arbitrary scale and v is a characteristic 
exponent with v > 0. With this selection, weights de- 
crease with increasing degree, in such a way that edges 
with the largest weight connect vertices with lowest de- 
gree. Obviously, an unweighted network will correspond 
to v = 0. 

(ii) Weighted BA tree II (WBAT-II): Weights as- 
signed according to their age in the network construction, 
namely 



l*-Jl 
N 



(4) 



where the node numbers i and j correspond to the par- 
ticular time step in which they were first introduced in 
the network. Since the degree of nodes decreases as 
k% oc {N/i) 1 / 2 during this process, where N is the size 
of the network, this selection with x > favors con- 
nection between unlike nodes and suppresses interactions 
between similar ones. 

The presence of weights in the network affects the dy- 
namics of the CP in the rate at which empty vertices be- 
come occupied. Thus, the rate at which an empty vertex 
i becomes occupied on contact with an occupied vertex 
j is now proportional to Xu%j. The proportionality of 
this rate with ujij implies that the creation of new parti- 
cles takes place with larger probability between vertices 
joined by a large weight edge. We can thus conclude that 
activity can in principle become trapped in isolated con- 
nected subsets of vertices, joined by large weight edges, 
subsets which are at the same time connected among 
them only through small weight edges, which transport 
activity only with very small probability. As we will see 
below, these correlated subsets of vertices will play the 
role on rare regions in the analysis of the CP in tree 
weighted networks. 



B. Heterogeneous mean-field analysis 

An analytical understanding of the CP (and in gen- 
eral of any dynamical process) on complex heterogeneous 
networks can be gained through the application of het- 
erogeneous mean-field (HMF) theory [111 13]. Thus ' in 
the case of symmetric weights depending on the degree at 
the edge's end points, i.e. cjy- = g(ki, kj), we can write, 
following Ref. [33|, the rate equation for the density of 
occupied vertices of degree fc, p k (t), namely 



Pk(t) 



E 



p k {i) + Xk[l - Pk (t)] 

g(k',k)P(k'\k) 



k'J2 Q g(k',q)P(q\k' 



■Pk>(t), (5) 



where P(k'\k) is the conditional probability that vertex k 
is connected to vertex k' [H| . The analysis of Eq. ([5]) for 
general functions g(k',k) and P(k'\k) presents notable 
mathematical difficulties. We will then focus in the case 
of uncorrelated networks, with P{k'\k) = k' P(k') / (k) 
and multiplicative weights g(k',k) = g(k)g(k') j33| . 



which correspond to the WBAT-I model, with g(k) = 
^/uJok~ v '. The former approximation is in this case es- 
sentially correct, since BA networks are very weakly cor- 
related [3!| . We are thus lead to the simplified rate equa- 
tion 



p k (t) = ~Pk(t) + X[l - p k (t)} 



k^pit) 



(6) 



where p(t) = ^2 k P(k)p k (t) is the total density of occu- 
pied nodes in the network. 

From Eq. ([S]), we see that p k = is always a solution. 
The conditions for the presence of non-zero steady states 
can be obtained by performing a linear stability analy- 
sis (40[. Neglecting higher order terms, Eq. © can be 
linearized as p k (t) ~ J2 k > L kk >p k '(t), with 



L kk ' 



Sk.k> + A 



k^Pik') 



(7) 



It is easy to see that the Jacobian matrix L kk i has a 
largest eigenvalue A = A— 1, associated to the eigenvector 
v k = k 1 ~ v . A nonzero steady state is only possible for 
A > 0, which translates in a critical threshold 



A,. 



1. 



(8) 



independent of the degree distribution and the weight 
pattern of the network 4 . 

More information on the active phase can be obtained 
by imposing the steady state condition p k (t) = 0, yield- 
ing the nonzero solution 



p k = 



Xk^p/ik 1 -") 
l + \k l -"p/{k 1 - 1 ') 



(9) 



Inserting this expression on the relation p = ^ fc P(k)p k , 
we can obtain a self-consistent equation for p in the 
continuous degree approximation, substituting sums by 
integrals and approximating the degree distribution by 
P(k) = 2m 2 k~ 3 with k G [m, oo] for BA networks. The 
form of this approximation depends on the value of v 
(see Appendix for a detailed calculation). So, for v = 0, 
which corresponds to a an unweighted BA network, we 
obtain a particle density, close to the critical point, given 
by 



A- 1 



hi 



' A-l ' 
. A - 



(10) 



We thus have a critical point A c = 1 and a mean-field 
exponent (3 = 1, with additional weak logarithmic cor- 
rections. For v > 0, on the other hand, we obtain, close 
to the critical point 



A-l, 



(11) 



4 We do not expect, however, to obtain this threshold in numerical 
simulations, since it is a non-universal parameter Illl. ll4)l . 
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recovering the homogeneous MF result exponent = 1. 

In order to find the decay of the order parameter at 
the critical point, we look at the time evolution of p(t), 
setting A — 1, namely 



(k) 



J2P(k)k Pk (t). (12) 



Performing a quasi-static approximation [15l |41| , substi- 
tuting pk(t) by the form given by Eq. §§§ with A = 1, 
and performing again a continuous degree approxima- 
tion, changing sums by integrals, we obtain a simple dif- 
ferential equation that can be easily solve. Thus, for the 
case v = we are led to (see Appendix) 



p(t) ~ [iln(i)] 



(13) 



which recovers the homogeneous mean-field result, p{i) ~ 
t~ a , with a = 1 plus additional logarithmic corrections. 
For v > 0, we obtain instead a full homogeneous mean- 
field result, a = 1. 

We conclude from this analysis that the CP on the 
weighted BA model should (at least at mean field level) 
show a critical behavior compatible with HMF exponent 
in the thermodynamic limit for v > 0, while weak loga- 
rithmic correction are expected for v = 0. 



C. Finite-size effects 

Recently in has been called attention on the finite-size 
effects in the CP on heterogeneous networks, which can 
drive the transition so strongly as to overcome the ther- 
modynamic limit behavior [U Q EH ■ It is easy to check 
that in the present case, these finite-size effects are very 
weak for v — or irrelevant for v > 0. In fact, following 
in a simplified form the reasoning in Ref. [15j, we can 
write the time evolution of the total density p(t) as 



p{t) = Y,Pk(t)p{k) 



Ap(t) - 



E 



fc2(l-*)p(fc) 



(14) 



where A = A — 1 and we have used Eq. ©. In the limit 
of very small particle density in a network of finite size 
N, when p -C (k 1 ^") /kl~ u , where k c is the maximum 
degree of cut-off in the network [42l ] , we can approximate 



Pit) - A P (t)-|^X>^>P(fc) 
= p(t) [A - X 2 p{t)g v ] , 
where we have defined the parameter 

(fc 2 ^-")) _ \ ln(iV) for v = 



const for v > 



(15) 



(16) 



The only nonzero solution of Eq. (fTS")) gives thus the finite 
size behavior 



A- 1 
A^7' 



(17) 



Form Eq. (|16p . we can see that g v is size independent 
for v > 0, while it depends weakly (logarithmically) on 
N for v = 0. We thus conclude that, even at the finite- 
size regime, the CP in weighted BA networks is fully 
described by the homogeneous mean-field critical point, 
with no strong corrections to scaling except for the pres- 
ence of weak, logarithmic ones at v = 0. 

Size effects also appear in the density decay with time 
at the critical point. We can estimate those effects by 
considering the time evolution of p(t) at A = 1 (A = 0) 
in a finite network. Again in the limit p < {k 1 ^) / 'k\~ v ', 
we can simply integrate Eq. (|15[) , to obtain the finite size 
density decay in the limit of large times as 



p(t) ~ [g v t\ 



(18) 



Again we observe that density decay in time is given ac- 
cording to the homogeneous mean-field theory for v > 0, 
with logarithmic size corrections at v = 0. 

To conclude this Section, we emphasize that the HMF 
solution presented here is expected to be correct only in 
the case of non-tree looped networks, since HMF theory is 
know to provide inaccurate predictions in tree topologies 
[32[ l35l . lia . |44| , and to be accurate on weighted networks 
only within certain limits [33l [45j . It is nevertheless in- 
structive to consider its results as a first order approxima- 
tion. Additionally, we must note that no HMF prediction 
is possible for the WBAT-II, due to the additive nature 
of the weight pattern. 



III. NUMERICAL SIMULATIONS 

In a practical numerical implementation of the CP, 
the neighbors of all nodes are stored in a dynamically 
growing table of vectors, in such a way that networks 
up to size N = 8 x 10 7 can be stored in 7.2 GB of com- 
puter memory. We follow the density decay of active sites 
by starting from fully occupied networks. Active nodes 
are selected randomly and deactivated with probability 
1/(1+A); otherwise in a randomly selected neighborhood 
a site activation occurs with probability [A/(l + A)]wjj. 
Time is updated by 1 Monte Carlo (MC) step after N a 
attempts, where N a denotes the number of active nodes 
in the previous time step. The simulations were run up 
to t = 2 22 - 2 26 MC time steps and have been averaged 
over 10 to 10 6 independent networks of different size. 
By the statistical sampling we used master-worker type 
of parallclization (MPI-C). This has been been repeated 
for several values of the control parameters: A, v and x. 
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FIG. 1: (Color online) Density decay (tp(t)) as a function 
of ln(£) for the CP on unweighted looped BA networks with 
m = 3 of size N = 8 x 10 7 . The different curves correspond 
to A = 1.2068, 1.24 (from bottom to top). Inset: Steady 
state density, showing agreement with HMF theory scaling. 
The full line shows a power-law fitting to the data points in 
the form -0.36(5)x°' 98(2) . 
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FIG. 2: (Color online) Density decay as a function of time for 
the CP on unweighted tree BA networks of size N — 8 x 10 7 . 
The different plots correspond to A = 2.756, 2.759 (from 
bottom to top). Right inset: The corresponding effective ex- 
ponents from a local slopes analysis (increasing A = from 
bottom to top). Left inset: Steady state density with power- 
law fitting of the form (|19|) for parameters in opposite order 
(from top to bottom) as in the main graph. 



A. Unweighted looped BA networks 

Wc fist consider the behavior of the CP on unweighted 
looped BA networks, generated with m = 3. Wc find a 
critical point around A c = 1.2068(1) and a critical ex- 
ponent f3 = 0.98(2) with logarithmic corrections, as ex- 
pected from Eq. (flU)) , see inset in Fig. [TJ The main plot 
of Fig. [1] presents the particle density decay for different 
values of A around the estimated critical point, rescaled 
according to the HMF prediction p(t) ~ l/[iln(i)]. This 
prediction is fairly well fulfilled by the simulations, as 
we can see from the almost straight linear part in the 
ln[tp(t)} vs. ln[ln(i)] curve for A = 1.2068, although fi- 
nite size corrections are noticeable. A linear fit in the 
region 4 < ln(ln(t)) < 7 for this value of A results in the 
slope —1.1(1), compatible with the HMF prediction. 

The behavior observed unweighted looped networks is 
therefore compatible with an standard critical point, rea- 
sonably described by HMF theory, showing a clear su- 
percritical phase reaching a stable plateau of activity for 
A > A c , an exponential activity decay for A < A c , and no 
apparent continuously varying power-law decays of activ- 
ity in the vicinity of the critical point, in agreement with 
the results reported in Ref. fl4l |. 



B. Unweighted BA trees 

We now investigate the effect of a loop-less tree struc- 
ture, again in the absence of weights. For small network 
sizes, N < 10 6 , numerical results do not yield to sim- 
ple interpretation. S-shapcd density decay curves can be 
seen on log-log plots, in such a way that the location of 
the transition point cannot be clearly determined. This 



observation suggests the presence of corrections to scal- 
ing that could be larger than those predicted in the sim- 
ple HMF analysis presented in Sec. IA 1 al In fact, as wc 
increase the network size, a phase transition with power- 
laws (PL) seems to emerge, but the corrections became 
negligible only for network sizes N > 10 7 nodes. In Fig. [2] 
we show the particle density decay on networks of size 
N = 8 x 10 7 in a narrow scaling region around the crit- 
ical point given by 2.756 < A < 2.759. We have done a 
careful analysis of the steady-state, with p leveling off to 
constant values, taking into account the effects of finite 
sizes. By assuming a critical transition (see left inset in 
Fig. [2]) in the form 



p{oo) = A\X - X c f , 



(19) 



we obtain a critical point A c = 2.757(1) with A = 0.287 
and ft = 1.20(5) by means of a least-squares fitting per- 
formed discarding the smallest values of A. Note that this 
exponent j3 ~ 1.2 is larger than the expected HMF result, 
an observation hinting towards a failure of the HMF the- 
ory for the CP on tree-like networks [32| • The main plot 
of Fig. [2] depicts the decay in time of the particle den- 
sity p(t) for different values of A around the estimated 
critical point. The curves shown exhibit a decay expo- 
nent apparently varying with the value of A at very large 
times. In order to explore in more detail the decay of the 
density functions, we have computed an effective decay 
exponents, defined as the local slope of p(t) as given by 



a e s(t) 



Hp(t)/ P (t')} 

ln(i/f) 



(20) 



where t and t' have been chosen in such a way that 
the discrete approximate of the derivative is sufficiently 
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FIG. 3: (Color online) Density decay as a function of time 
for the CP on weighted BA trees generated with the WBAT- 
I model with exponent v — 1.5. Network size N — 10 5 . 
Different curves correspond to A = 160, 156, 154, 149, 148, 
147, 146, 145, 144.7, 144.2, 144, 143.5, 143, 142, 140 (from 
top to bottom). Inset: Steady state density. 



FIG. 4: (Color online) Density decay as a function of time p(t) 
(lowest branch curves) and survival probability P(t) (highest 
branch curves) for seed simulations of the CP on weighted 
BA trees generated with the WBAT-I model with exponent 
v = 1.5. Curves correspond to A = 145, 144, 143, 140 (top to 
bottom for each branch). Network size N = 10°. 



smooth, see right inset in Fig. [5J From this figure, we 
observe that the decay exponents are indeed different, 
smaller than the HMF prediction, which would lead to an 
effective exponent larger than one in a log-log plot. More- 
over, the exponents thus computed are not well defined 
for all values of A, since the plots of a e g(t) do not show 
in general well-defined plateaus except for very small re- 
gions at very large times. 

The conclusion we draw from the analysis of non- 
weighted BA trees is that these network structures do not 
show evidence for the presence of explicit RR effects, but 
are compatible with a standard non-equilibrium phase 
transition. Evidences for RR effects, and in particular 
for a smeared transition, are however compelling in the 
case of weighted BA trees, as we will see below. 



C. Weighted BA trees with multiplicative weights: 
WBAT-I model 

We next investigate the possibility of inducing RR ef- 
fects in the activity propagation of the CP by adding 
to the BA tree networks a weight structure wich will act 
against the "stronger gets stronger" behavior of B A mod- 
els by suppressing the dominancy of strongly connected 
nodes. 

In the first place, we consider the multiplicative weight- 
ing strategy defined by model WBAT-I. In this case, we 
have as a theoretical guideline for the behavior of the CP 
the results from the HMF analysis presented in Sec. Ill BL 
which suggest a fully homogeneous mean-field behavior, 
with j3 = a = 1. We note that departures from this 
prediction are however not unexpected, given the tree- 
like and weighted structure of the network model under 
consideration. 



For a small weight exponent v = 1/2 we find a phase 
transition at A c = 8.62(1), characterized by the mean- 
field type of density decay: p(t) cx 1/t (not shown) 
[4fj| . For larger values of v, however, generic, continu- 
osly changing power-laws can be observed. In particular 
for v = 1.5 simulations performed on networks of size 
N = 10 5 suggest the presence of a GP as shown on Fig. [3] 
Considering the steady state value of the density, we can 
observe a very smooth approach to zero. A power-law 
fitting to the form p = A\ \ — X c \^ results is A c = 143(1) 
and f3 = 1.3(1), as shown in the inset of Fig. [3] The value 
of /3 obtained in this fitting is again in disagreement with 
the HMF prediction from Sec. Ill Bl 

We have also performed spreading simulations, start- 
ing from a very small initial seed of active sites @ . These 
simulations were performed on networks with N = 10° 
nodes up to t = 4 X 10 5 MC time steps. Besides the 
density pit), we measured the survival probability P(t) 
[1], averaged over 10 8 independent runs on 10 3 differ- 
ent network configurations. As Fig. [4] shows, not only 
the density, but also the survival probability P(t) curves 
exhibit slow, power-like behavior for times larger than 
t > 2000. In the initial time region t < 1000, both P(t) 
and p(t) decay exponentially. 

Simulations performed for larger network sizes suggest, 
however, that the apparent PLs observed in fact disap- 
pear in the thermodynamic limit of large network size. 
In Fig. Owe plot the particle density decay as a function 
of time for A = 140 and 144 corresponding to increasing 
network sizes from N = 10 5 up to N = 10 6 . As can be 
seen from this plot, the apparent power laws, fitted by 
dashed lines, that can be observed for N < 10 5 , turn into 
saturation plateaus for N = 10 6 at very late times. 

This observation suggests that what we are actually 
observing in the WBAT-I model is the tail behavior of 
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FIG. 5: (Color online) Density decay as a function of time p(t) 
for the CP on weighted BA trees with a multiplicative weight- 
ing scheme (WBAT-I) with exponent v = 1.5. Plots corre- 
spond to two sets of A (upper branch: A = 144, lower branch 
A = 140) at different network sizes N. Dashed lines represent 
PL fittings. Inset: Initial time region of the same data, show- 
ing an stretched exponential behavior near the transition. 



a smeared transition [2l| . where the density has a con- 
stant nonzero value p s t = p(oo) in the long time limit 
and p(t) approaches this in a PL manner, see Eq. ([TJ. 
Additional evidence of this picture come from the analy- 
sis of the initial time decay of the particle density, with 
can be fitted to an stretched exponential of the form 
p oc exp(t -0 - 91 ) (see inset on Fig. [5]). This behavior is 
very close to an exponential decay. However, comparing 
with the CP with correlated disorder on regular lattices, 
at intermediate times of the smeared transition we should 
observe a stretched exponential of the form in Eq. ([2]) . A 
value d r 3> z in Eq. <j2j) would fit to the observed behav- 
ior, an expectation that is quite reasonable in a random 
network with infinite topological dimensions. 

The conclusion extracted from the analysis of this 
model is that the presence of a smeared phase transition 
is quite probable. A smeared transition can explain the 
saturation in the thermodynamic N — > oo limit, due to 
the presence of correlated subspaces of connected vertices 
of small degree, and correspondingly joined by relatively 
large weights. These subspaces play the role of RRs, that 
can remain active, leading to an initial stretched expo- 
ncntia l dy namics crossing over to power-law for longer 
times 
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In order to visualize the presence of these RRs, we have 
performed a percolation analysis psj - IHoj of our WBAT- 
I networks. We consider a network of a given size N, 
and delete all the edges with a weight smaller than a 
threshold wth- For small values of w t h, many edges re- 
main in the system, and they form a connected network 
with a single cluster encompassing almost all the vertices 
in the network. When increasing the value of Wth, the 
network breaks down into smaller subnetworks of con- 
nected edges, joined by weights larger than w t h- The 




FIG. 6: (Color online) Size Si of the 5 largest clusters in a 
percolation analysis of the WBAT-I model with v = 1.5 for 
wth = lOOojmin (hollow symbols) and Wth = 1000oj m i n (full 
symbols), where oj m i n is the minimum weight in the network. 
The size of all components grows linearly with network size 
N, and is therefore infinite in the thermodynamic limit. 



expectation in the latter case from standard percolation 
in networks is to observe the presence of a largest cluster 
(the giant component) with a size scaling with N, while 
the rest of the components should have a size scaling log- 
In Fig. [6] we plot the 
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arithmically with N at most 
average size (measured as the number of vertices) of the 
5 largest components corresponding to a given threshold 
Wth- For the values of Wth considered, the networks break 
down in a number of components. However, the size of 
the largest ones grows linearly with the network size N, 
at odds with the expectations from a standard percola- 
tion transition. These clusters, which can become arbi- 
trarily large in the thermodynamic limit, play the role 
of correlated RRs, sustaining independently activity and 
smearing down the phase transition. 

This sort of analys allows also to explain why looped 
weighted networks do not exhibit RR effects. In this case 
(data not shown), the largest component scales again lin- 
early with the network size, but the next largest compo- 
nents scale in size logarithmically. They thus become ir- 
relevant in the thermodynamic limit and cannot sustain 
activity at very large times. 



D. Weighted BA trees with age-dependent weights: 
WBAT-II model 



We finally consider the behavior of the CP on weighted 
tree BA networks constructed with the age-dependent 
weighting scheme, defined by the model WBAT-II. In this 
case, the weights are anti-correlatcd with the degree, in 
opposition with the previous WBAT-I, in which large de- 
gree vertices had small weight, with small degree vertices 
had large weight. Here instead, weight is larger in edges 
connecting vertices with different degrees, while edges 
between vertices of similar degree have a small weight. 

Selecting an exponent for the weight strength x = 2 
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FIG. 7: (Color online) Density decay as a function of time 
p(t) for the CP on weighted BA trees with a age-dependent 
weighting scheme (WBAT-II) with exponent x = 2. Network 
size TV = 10 5 . Different curves correspond to A = 6.75, 6.8, 
6.85, 6.87, 6.9, 6.92, 6.94, 6.96, 6.98, 7, 7.04, 7.1, 7,2, 7.4, 8.5, 
9, 12, 15 (from top to bottom). Inset: Corresponding local 
slopes for A = 6.9, 6.92, 6.94 (from bottom to top). 
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FIG. 8: (Color online) Finite-size scaling analysis of the den- 
sity decay exponent for A = 6.75 (triangles), A = 6.8 (boxes), 
A = 6.82 (triangles), A = 6.85 (bullets), A = 6.9 (rhombes) in 
the CP on weighted BA trees with a age-dependent weighting 
scheme (WBAT-II) with exponent x — 2. Top inset: p(t) for 
A = 6.82 (N = 10 6 , N = 4 x 10 5 , TV = 10 5 from top to bot- 
tom). Bottom inset: Initial time density near the transition. 



and not very large network sizes, N = 10 , generic power- 
laws can be observed for a large range of A of values, as 
shown on Fig. [71 which exhibit a very strong similarity 
with those found in the multiplicative WBAT-I model, 
see Fig. G3 An analysis of the steady state density at 
large times in the region 6.8 < A < 9 can be fitted again 
to the canonical form p = A\X — \ c \@ , with (3 = 1.35(5), 
A = 0.02(5) and A c = 6.80(2). These values are again 
in discrepancy with HMF expectation, but this fact is 
altogether not surprising, since we do not have a proper 
HMF solution for this model, and we can only proceed 
to compare with the homogeneous result observed for 
the WBAT-I model. The inset in Fig. [7J shows the ef- 
fective exponents a e s(t) as computed using the formula 
in Eq. (|20[) . From the figure we observe the presence of 
reasonably large plateaus that confirm the presence of 
well-defined PLs with exponent varying with the value of 
A. 

Simulations over larger networks hint again towards 
the disappearance of the PLs observed in smaller network 
sizes (see top inset of Fig. [5]). In FigJS] we plot the value 
of the decay exponent for larger systems by performing 
a power-law regression in the last two decades of time. 
This seems to indicate that as a function of the network 
size for A > 6.8, the estimated exponent apparently tends 
to zero, suggesting the presence of an active phase in the 
thermodynamic limit. On the other hand, for A < 6.8 the 
decay exponent does not change too much with N and 
seems to stay a > 1 up to the largest sizes N = 4 x 10 6 
we could reach by extensive simulations. 

This kind of behavior can again be understood in terms 
of a smeared phase transition. Further evidence of this 
fact comes from the analysis of the initial decay of the 
activity density, which again can be fitted to a stretched 




FIG. 9: (Color online) Size Si of the 5 largest clusters in a 
percolation analysis of the WBAT-II model with x — 2 for 
Uth = lOOcJmin (hollow symbols) and wth = 1000oj m i n (full 
symbols), where u m i n is the minimum weight in the network. 
Curves in this last case have been shifted downwards for the 
sake of clarity. 



exponential of the form p cx cxp(i~°- 87 ( 5 )) (see bottom 
inset of FigjS]), in agreement with optimal fluctuation 
theory in an infinite dimensional network, Eq. (|2[). In 
this case, the corresponding correlated subspacc will be 
formed by a subset of connected vertices with alternating 
small and large degree, and consequently joined by large 
weight edges, which would remain active in the thermo- 
dynamic limit giving rise to the characteristic RR effects 
observed in simulations. The nature of the RRs in this 
model can also be visualized by means of a percolation 
analysis. In Fig. Owe plot the average size of the 5 largest 
components in percolation experiments for different val- 
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FIG. 10: (Color online) Density decay as a function of time 
p(t) for the CP on weighted BA trees with a age-dependent 
weighting scheme (WBAT-II) with exponent x = 3. Different 
curves correspond to A = 9.4, 9.5, 9.6, 9.65, 9.7, 9.8, 9.85, 
(from bottom to top). Left inset: The corresponding effective 
exponents from a local slopes analysis for A = 9.6, 9.65, 9.7, 
9.8, 9.85 (from bottom to top). Right inset: initial, stretched 
exponential behavior. 



ues of the weight threshold Wth- Again, all five scale 
linearly with network size, confirming their role as corre- 
lated RRs in the thermodynamic limit. 

This sort of behavior is not exclusive of the value x = 
2 considered in the simulations above, but can also be 
observed for other values of x, see Fig. [TO] The initial 
time dependence for t < 30 MCS follows in this case 
p oc exp(*-°- 87 ( 5 )) (see right inset of Fig. [TO}. 



IV. CONCLUSIONS 

The heterogeneous and disordered pattern of connec- 
tions in a complex network can induce considerable and 
surprising differences in dynamical processes running on 
top of them, as compared with their outcome on or- 
dered or homogeneous systems. Building on recent re- 
sults concerning the additional effects of intrinsic ran- 
dom, quenched disorder superimposed over a network 
structure, in this paper we have shown that topology- 
induced disorder can induce rare-region effects on the 
contact process. In particular, we have investigated the 
effects of a tree topology, in which no global loops are 
present in the network (there is a unique path between 
any two vertices) and of a correlated weight pattern as- 
signed to the set of edges. Both elements have been 
show to induce a slowing down in dynamics as simple 
as diffusion. By means of large scale numerical simula- 
tions on tree-like weighted networks generated using the 
Barabasi- Albert model, we have shown that a tree-like 
topology is not by itself enough to induce a slow dynam- 
ics or rare-region effects. However, the combined effects 
of a tree structure and a correlated weight pattern are in- 



deed capable of triggering a very noticeable slow dynam- 
ics, which is compatible with the presence of rare-region 
effects. The two different weighting schemes we have 
considered for small system sizes suggested the appear- 
ance of Griffiths phases, characterized by a power-law dy- 
namics of the activity and the survival probability above 
the clean critical point, with a varying decay exponents. 
More accurate simulations performed on larger network 
sizes indicate that, instead of true Griffiths phases, the 
behavior of the contact process on tree weighted networks 
is better explained in terms of a fully smeared transi- 
tion, characterized by an initial time decay of the activity 
density with the form of an stretched exponential, whose 
presence can be understood by optimal fluctuation theory 
arguments applied on a network with infinite topological 
dimension. The smearing of the transition can be argued 
to be due to the presence of a correlated subspace of con- 
nected vertices, which are joined by large weight edges, 
playing the role of rare-regions and trapping activity at 
different large time scales. 

To sum up, we have shown here that topological dis- 
order is capable to induce rare-region effects on the con- 
tact process on complex networks. In the particular case 
we have considered (weighted tree networks), those rare- 
region effects translate in the smearing of the critical 
phase transition. Our work opens the path to a deeper 
investigation of this issue, focusing in particular on dis- 
covering those new topological ingredients which could 
induce fully-Hedged Griffiths phase originating on topo- 
logical grounds. 
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Appendix A: Analytic solution of the HFM 
equations for WBAT-I model 

In this section we develop the full HMF solution of 
the WBAT-I model on uncorrelated networks, sketched 
in Section iHBl 
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1. Steady state density 

The starting point is Eq. ([6]), namely 

,k l - v p{t) 



Pk(t) = -p fc (t) + A[l -p k {t)]- 



(Al) 



which, in the steady state regime, leads to the non-zero 
solution for the partial density of particles in vertices of 
degree k 



Pk 



Xk 1 -"p/(k 1 - v ) 
1 + Xk 1 -" p/{k l - v ) 



(A2) 



From this expression, p, the total density of particles, can 
be computed self-consistently, noticing that 



E p ( fc ) 



Pk 



k 

2m 2 



k 

Xk~ 2 -" 'p/ik 1 -" 



1 + X^-vp/ik 1 -"} 



1 + Xk^p/ik 1 -") 

dk : (A3) 



where in the last expression we have used the degree dis- 
tribution in its continuous degree form, P{k) = 2m 2 fc~ 3 , 
and substituted sums by integrals. Performing the inte- 
gral in Eq. (|A3|) leads to an equation that can be solved 
to obtain p as a function of A. The form of this equation 
depends on the value of v. Let us consider separately the 
different cases. 



where F[a, 6, c, z] is the Gauss hypergeometric function 
[521 ] . In order to evaluate the critical behavior of the den- 
sity in the thermodynamic (infinite network size) lime, 
we can invert the previous expression using asymptotic 
expansion of the Gauss hypergeometric function [52| . 
namely 



F 



1, a, 1 + a, ■ 



1 



T(q + l)r (1 - a) + a^(-l) 



(A8) 



n=l 



For small p (small z), since 2/(1 — v) > 2 for v > 0, the 
most relevant terms are z and z 2 . Thus, at leading order, 
the have 



1 / Xpm 



l-v 



2 Xpm 1 



(A9) 



Using now the fact that, in the continuous degree ap- 
proximation, = 2m 1 ~' y /(l + v), we are led to the 
expression 



(A-l), 



' (l + i/) 2 A 2 
recovering the asymptotic expression Eq. (fTTj) . 



(A10) 



a. v = 



In this case, Eq. (|A3|) takes dc form 
"°° Xk- 2 p/(k) 



p = 2m 
= Xp 



l + Xkp/(k) 



dk 



2 \x P 



(A4) 
(A5) 



where in the last expression we have used (k) = 2m. 
Eq. (|A5[) provides an implicit equation for p which, in 
the limit p — > 0, close to the critical point, leads to the 
explicit solution 



2 A-l 



(A6) 



which corresponds to the asymptotic behavior in 
Eq. (HDD- 



b. < v < 1 

In this case, performing the integral for k in Eq. (|A3I) . 
we obtain the self-consistent equation 



p = F 



-, 1 



l-v' l-v' \ (k 1 -") 



Xpm 1 



(A7) 



C. V > 1 



Now, the self-consistent equation for p is 



p = XpT 



2v 



F 



v + 1 v + 1 Xm x ~ v p 

E— ^,1 1 



(All) 

where F[a, &, c, z] is the regularized Gauss hypergeomet- 
ric function. Using the series expansion for small z [52| 

T(a + l)F[l,a,a+ 1,-z] = [ ~ ' * , (A12) 



n=0 



n + a 



the self-consistent equation at leading order in p takes 
the form 



p~ Xp 



v — 1 Am 
1 — - , ... p 



from where it follows that 



+ 1 (fc 1 "") 



r(A-i), 



(^-1)A 2 ' 

from which the scaling form Eq. pip ensues. 



(A13) 



(A14) 
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2. Density decay at criticality 

The decay of the particle density at criticality can be 
obtained from Eq. (fT2")l . namely 



p( t ) = j2p(k)m = -jtlY, p w k p^ 

k k 

, s f°° Xk- 1 -" pit) I Ik 1 -") „ , k , 

where in the last expression we have performed a quasi- 
stationary approximation, substituting pk (t) by the func- 
tional form in Eq. (|A10[) . and passed to the continuous 
degree approximation. 



a. v = 



Eq. (|A15[) takes the form 

o(*) = ^U i + AMt)/w rffc 
- >('>Vi' 2 



2 -"V' ' p(t)' 

This equation can be integrated, yielding 



(A16) 



(A17) 



where li(z) is the logarithmic integral function [52j. Using 
the asymptotic expansion li(-z) ~ zj ln(z) , we obtain the 
density decay in the large t limit limit given by Eq. (fTc 



b. v>Q 



In the case < v < 1, integration of Eq. (|A15|) leads 



to 



p{t) ~ 



l + v 
l — i/ 



1 + 1/ f\pm 1 -" x 



i - 1/' v ( fcl_,y > 

(A18) 

Expanding for small density, using Eq. (|A8[) , and noticing 
that at (v + 1)/(1 — v) > 1, we have at leading order in 
P, 



m 



l + i/ Am 1- " , 



(A19) 



leading to the density decay, in the infinite network size 
limit, p(t) ~ t- 1 . 

An analogous calculation, implying the regularized 
Gauss hypergeometric functions, can be performed for 
v > 1, leading to the same mean- field behavior for the 
density decay. 
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